function f = expweightedavg(DATA , nu , drop)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function calculates the exponentially-weighted average of all past
% data, i.e.
%    y_t = \nu\sum_{i=0}^{t-1}(1-nu)^i x_{t-i}
% for each t.
%
% INPUT:
%     DATA - T by K matrix
%     nu   - weights on most recent observation
%     drop - burin-in period [OPTIONAL]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%--------------------------------------------------------------------------
% Check inputs
if nargin < 2
    error('2 inputs required!')
end

switch nargin
    case 2
        drop = 0;
end

%--------------------------------------------------------------------------
% Get size
[T , K] = size(DATA);

%--------------------------------------------------------------------------
% To speed up the calculation of exponential-weighted average, we use a
% cumsum approach with adjustment for exponential term
% However, for very long series, the first observation in the algorithm
% will be reduced to zero in Matlab precision, which results a 0/0
% situation and makes the algorithm very inaccurate. We adjust for this
% problem by multiplying a reasonable exponential term to make it not too
% small. The following section determines how much an exponential term
% should be added based on Matlab precision
nmin = realmin('double');
nmax = realmax('double');

data_min = min(min(abs(DATA(1 : floor(T/2) , :))));
data_max = max(max(abs(DATA(floor(T/2) + 1 : end , :))));

if data_min < nmin
    data_min = nmin;
end
% This is the minimal T_0 needed to be added to exponential
T_0_lb = max(ceil(T - 1 - log(nmin/(data_min * nu))/log(1 - nu)) , 0);

% This is the maximal T_0 needed to be added to exponential
T_0_ub = min(-log(nmax/(data_max * nu))/log(1 - nu) , T);

% The use of T_0 is to prevent numerical error in Matlab when dealing with too small a number
T_0 = round((T_0_lb + T_0_ub)/2); 

%--------------------------------------------------------------------------

% Calculate the exponentially-weighted average of past data
f = cumsum( nu * DATA .* repmat((1 - nu).^(T - 1 - T_0 : -1 : 0 - T_0)' , 1 , K) )./repmat((1 - nu).^(T - 1 - T_0: -1 : 0 - T_0)' , 1 , K);

% Drop burn-in periods if needed
f = f(drop + 1 : T , :);
end